function  [D_CS_avg, D_Ext_avg, D_Mis1_avg, D_Mis2_avg, D_Mis_avg, D_SW_avg, m_theta, m_D_E, m_E, sd_theta, sd_D_E, sd_E, corr_theta_D_E, corr_theta_E, D_Mis_SumStats, D_Mis_SumStats_1, D_Mis_SumStats_2, D_Mis_SumStats_3, D_Mis_SumStats_4] = Total_Welfare_decomp_WTP_EEgap_parfor(weight,param,param_FE,param_FEdemo,id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,MT,price, price_denom, rebate_estar_denom, estar_denom, kwh, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom)
                                                                                                                                                                                                                                                    
    param_FE_k=param_FE;
	param_FEdemo_k=param_FEdemo;

dim_k=length(weight)
%parfor k=1:dim_k
for k=1:dim_k
    
	param_k=param{k};
    eta_k(k)=param_k(1);
    theta_k(k)=param_k(3);

%The social price of electricity is computed as follows: p_elec_ext = elec+p_elec+t_pigou;

%Welfare evaluated at the baseline 
d_ext_J = repmat(d_ext',J,1);
t_pigou_0=t_pigou*0;
[CS_0_i{k}, CS_0_avg(k), E_0_i{k}, E_0_avg(k),  W_0_i{k}, W_0_avg(k), Error_0_i{k}, Error_0_avg(k), ext_0_i{k}, ext_0_avg(k), gvt_0_i{k}, gvt_0_avg(k), SW_0_i{k}, SW_0_avg(k), P_0_decision{k}, P_0_decision_j{k}] = ...
Welfare_decomp_k_WTP_EEgap_parfor(param_k,param_FE_k,param_FEdemo_k,id_ij,J,p_elec,t_pigou_0,d_ext_J,elec,lifetime,r,MCPF,MT,price, price_denom, rebate_estar_denom, estar_denom, kwh, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom);

%Welfare evaluated at the social price of electricity 
[CS_i{k}, CS_avg(k), E_i{k}, E_avg(k),  W_i{k}, W_avg(k), Error_i{k}, Error_avg(k), ext_i{k}, ext_avg(k), gvt_i{k}, gvt_avg(k), SW_i{k}, SW_avg(k), P_decision{k}, P_decision_j{k}] = ...
Welfare_decomp_k_WTP_EEgap_parfor(param_k,param_FE_k,param_FEdemo_k,id_ij,J,p_elec,t_pigou,d_ext_J,elec,lifetime,r,MCPF,MT,price, price_denom, rebate_estar_denom, estar_denom, kwh, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom);


D_CS_ik(:,k)   = CS_i{k}-CS_0_i{k};
D_Ext_ik(:,k)  = d_ext.*(E_i{k}-E_0_i{k});
D_Mis1_ik(:,k) = (theta_k(k)-1)*elec.*(E_i{k}-E_0_i{k});
D_Mis2_ik(:,k) = theta_k(k)*t_pigou.*E_i{k};
D_Mis_ik(:,k)  = (theta_k(k)-1)*elec.*(E_i{k}-E_0_i{k}) + theta_k(k)*t_pigou.*E_i{k};
D_SW_ik(:,k)   = SW_i{k}-SW_0_i{k};

%We don't consider heterogeneity at the individual (i) level with respect
%to the externality (d_ext), electricity prices (elec), and t_pigou
%emission factors (t_pigou)
D_CS_avgk(k)   = CS_avg(k)-CS_0_avg(k);
D_Ext_avgk(k)  = mean(mean(d_ext))*(E_avg(k)-E_0_avg(k));
D_Mis1_avgk(k) = (theta_k(k)-1)*mean(elec)*(E_avg(k)-E_0_avg(k));
D_Mis2_avgk(k) = theta_k(k)*t_pigou*E_avg(k);
D_Mis_avgk(k)  = (theta_k(k)-1)*mean(elec)*(E_avg(k)-E_0_avg(k)) + theta_k(k)*t_pigou*E_avg(k);
D_SW_avgk(k)   = SW_avg(k)-SW_0_avg(k);

%Check_D_SW_avgk(k) = D_CS_avgk(k) + D_Mis_avgk(k) - D_Ext_avgk(k);
%Check_SW_avgk(k) = CS_avg(k) + Error_avg(k) - ext_avg(k) + gvt_avg(k);
%Check_SW_0_avgk(k) = CS_0_avg(k) + Error_0_avg(k) - ext_0_avg(k) + gvt_0_avg(k);
%Check_Mis_avg(k) = Error_avg(k)+ gvt_avg(k)
%Check_Mis_0_avg(k) = Error_0_avg(k)+ gvt_0_avg(k)

end

    D_CS   = sum(weight'.*sum(D_CS_ik,1)/MT);
    D_Ext  = sum(weight'.*sum(D_Ext_ik,1)/MT);
    D_Mis1 = sum(weight'.*sum(D_Mis1_ik,1)/MT);
    D_Mis2 = sum(weight'.*sum(D_Mis2_ik,1)/MT);
    D_Mis  = sum(weight'.*sum(D_Mis_ik,1)/MT);
	D_SW   = sum(weight'.*sum(D_SW_ik,1)/MT);

    D_CS_avg  = sum(weight'.*D_CS_avgk);
    D_Ext_avg = sum(weight'.*D_Ext_avgk);
    D_Mis1_avg = sum(weight'.*D_Mis1_avgk);
    D_Mis2_avg = sum(weight'.*D_Mis2_avgk);
    D_Mis_avg = sum(weight'.*D_Mis_avgk);
	D_SW_avg  = sum(weight'.*D_SW_avgk);

%Compute sufficient stats
m_theta        = sum(weight'.*theta_k); 
m_D_E          = sum(weight'.*(E_avg-E_0_avg));
m_E            = sum(weight'.*E_avg);
sd_theta       = (sum(weight'.*(theta_k-m_theta).^2))^0.5;
sd_D_E         = (sum(weight'.*((E_avg-E_0_avg)-m_D_E).^2))^0.5;
sd_E           = (sum(weight'.*(E_avg-m_E).^2))^0.5;
cov_theta_D_E = sum(weight'.*((theta_k-m_theta).*((E_avg-E_0_avg)-m_D_E)));
cov_theta_E   = sum(weight'.*((theta_k-m_theta).*(E_avg-m_E)));
corr_theta_D_E = cov_theta_D_E/(sd_theta*sd_D_E) ;
corr_theta_E   = cov_theta_E/(sd_theta*sd_E);

%Reconstruct the welfare component based on sufficient statistics
D_Mis_SumStats_1 = mean(elec)*(m_theta-1)*m_D_E;
D_Mis_SumStats_2 = mean(elec)*corr_theta_D_E*sd_D_E*sd_theta;
D_Mis_SumStats_3 = t_pigou*m_theta*m_E;
D_Mis_SumStats_4 = t_pigou*corr_theta_E*sd_theta*sd_E;

%Compare whether the sufficient stats approach matches: D_Mis_SumStats vs D_Mis
D_Mis_SumStats = D_Mis_SumStats_1 + D_Mis_SumStats_2 + D_Mis_SumStats_3 + D_Mis_SumStats_4;

%Compare the sufficient stats approach for the heterogeneity components
%only:
D_Mis_SumStats_2_check = D_Mis1_avg-D_Mis_SumStats_1;
D_Mis_SumStats_4_check = D_Mis2_avg-D_Mis_SumStats_3;

end
